#ifndef VECTOR_H_
#define VECTOR_H_

#include <ostream>
#include <cmath>

using namespace std;

namespace geom {

template<class T>
class Vector4 {
public:
	T x;
	T y;
	T z;
	T w;
	Vector4(T x = 0, T y = 0, T z = 0, T w = 1) :
			x(x), y(y), z(z), w(w) {
	}

	template<class W>
	Vector4(const Vector4<W> &v) :
			x(v.x), y(v.y), z(v.z), w(v.w) {
	}

	T* getArray() {
		return &x;
	}

	const T* getArray() const {
		return &x;
	}

	T lengthSquare() const {
		return x * x + y * y + z * z;
	}

	T length() const {
		return sqrt(lengthSquare());
	}

	static T distance(Vector4<T> a, Vector4<T> b) {
		T dx = a.x - b.x;
		T dy = a.y - b.y;
		return sqrt(dx * dx + dy * dy);
	}

	void normalize() {
		T len = length();
		x /= len;
		y /= len;
		z /= len;
	}

	Vector4<T>& operator +=(Vector4<T> v) {
		x += v.x;
		y += v.y;
		z += v.z;
		return *this;
	}

	Vector4<T>& operator -=(Vector4<T> v) {
		x -= v.x;
		y -= v.y;
		z -= v.z;
		return *this;
	}

	Vector4<T> operator +(Vector4<T> v) const {
		return Vector4f(x + v.x, y + v.y, z + v.z, w);
	}

	Vector4<T> operator -(Vector4<T> v) const {
		return Vector4<T>(x - v.x, y - v.y, z - v.z, w);
	}

	Vector4<T> operator -() const {
		return Vector4<T>(-x, -y, -z, w);
	}

	float operator ^(Vector4<T> v) const {
		return x * v.x + y * v.y + z * v.z;
	}

	Vector4<T> cross(Vector4<T> v) const {
		return Vector4<T>(y * v.z - z * v.y, z * v.x - x * v.z,
				x * v.y - y * v.x);
	}

	Vector4<T> operator*(T v) const {
		return Vector4<T>(x * v, y * v, z * v, w);
	}

	friend Vector4<T> operator*(T s, Vector4<T> v) {
		return Vector4<T>(v.x * s, v.y * s, v.z * s, v.w);
	}

	friend ostream& operator<<(ostream& os, Vector4<T> v) {
		os << "(" << v.x << ", " << v.y << ", " << v.z << ")";
		return os;
	}

	static const Vector4<T> ZEROS, X_AXIS, Y_AXIS, Z_AXIS;
};

typedef Vector4<float> Vector4f;
typedef Vector4<double> Vector4d;

} /* namespace geom */

#endif /* VECTOR_H_ */
